Propagation of a wave through a medium can be described by a “wave equation.” The wave equation can typically be solved to determine the magnitude or other property of the wave as a function of time and space.
Geophysical surveys may be performed to explore the subsurface of the Earth. For example, seismic or electromagnetic surveys may be performed where an acoustic or electromagnetic (EM) wave is transmitted from a source into the subsurface and detectors are used to detect a response. Based on the detected response, properties of the earth may be determined. For example, the response may comprise time-series data with features associated with interfaces in the subsurface.
In seismic reflection surveying, high amplitude events in seismic data may be related to interfaces across which rock properties change abruptly and reflect a significant amount of energy back to the detectors. Therefore the high amplitudes in the data may be indicative of subsurface interfaces.
In typical processing of seismic data, the time between the transmission of the seismic wave and the detection of the high-amplitude event is often regarded to be the travel time of a seismic wave to the interface and back due to a reflection at the midpoint between source and receiver. Travel time is therefore typically considered a proxy for depth, and plotting of seismic data as a function of travel time may therefore indicate the depth relationship of amplitude features and hence subsurface interfaces.
In reality however, the seismic velocity of the subsurface layers through which the wave has passed is not constant (typically increases with depth), and therefore in order to obtain a more realistic impression of the location of interfaces, the data are depth converted, i.e. plotted as a function of depth, based on a pre-determined velocity profile or velocity model.
Data for source-receiver pairs are plotted at a mid-point therebetween. The data from all such pairs can therefore appear in a time or depth record at that midpoint, with amplitude events appearing as though they derive from a normal incidence source-receiver arrangement. Whilst this can give a suitable representation of the Earth's structure where layers are horizontal and planar, it will be appreciated that a seismic wave may be refracted along its propagation path, be reflected from sloping interfaces and interact with different layers in the incident and reflected directions. Therefore, the travel time data and depth converted data records may suggest a false geometry. For example, a high amplitude event may occur in the data at a particular depth position where there is in reality no reflector at that depth position.
This is a known issue in seismic data processing and in the interpretation of seismic reflection data, particularly in geologically complex structures. Correction of the data may be attempted by a process of migration. Seismic migration serves to reposition seismic data to their representative, geometrically faithful, positions in the data set. Migration algorithms have been developed for migrating seismic data. It is useful to understand the wave propagation of the seismic wave through the subsurface in order to perform migration. In simple terms, by tracing the route that the wave has taken from the source to the reflector and back to the surface, it is possible to identify the receivers at which the energy associated with that reflector has arrived and the time of arrival of that energy. In order to do this, a velocity model of the subsurface is required. This may be obtained from a model building package and may be based on well data or obtained in other ways, for example estimated, perhaps based on regional geological knowledge.
A further complicating factor is anisotropy of the subsurface. For example, particular rock types may have a strong directionality of propagation. The seismic velocity in different (e.g. orthogonal) directions may be significantly different, and may differ between layers.
Full wave-field migration techniques have been developed to model the development of the wave in time and at closely spaced positions spatially in 3D across a region of the subsurface. Such techniques may involve solving a wave equation of the kind mentioned above, to obtain the wave field at the different locations, with respect to time relative to a source signal emission at different locations.
These migration algorithms have as an input standard seismic data. The seismic data could be organized as common shot gathers, or in some other domain, such as plane wave or delayed shot domains. The seismic data may be provided in a transformed domain (compared to the originally acquired seismic data), for example a combined domain of linear phase for shot location for delayed shot migration, in which the combination of the data is performed by tau-p transform to the shot location. The tau-p transform may also be performed for receiver location to perform a plane wave migration. The wavefield solution of the wave equation provides a set of complex values at each target location/point, providing for example amplitude, phase, wave number, etc. These value are a crucial input components for the migration algorithms.
Many different wave equation techniques have been developed, and it is of interest to develop a wave equation technique for use in seismic migration in the case of an anisotropic subsurface. Some known migration techniques are now described in more detail.
Nowadays, accurate seismic anisotropic models are recognized as key points in imaging deep challenge areas for complex structures, because they describe the discrepancy of directional propagation speed of seismic waves, which is much more realistic for overburden geological structure such as bedding sediments with fractures. The tilted transverse isotropic (TTI) or tilted orthorhombic anisotropy (TOA) are generally necessary for imaging such complicated anisotropic structures (Zhang and Zhang, 2008; Zhang and Zhang, 2011; Fowler and King, 2011). Given an accurate anisotropic model, migration algorithms are destined to provide the artifact-free subsurface images. This may require applications of least-squares migration algorithms (Lambaré et al., 1992; Nemeth et al., 1999) and it may become even more complicated when working on imaging the elastic model (Jin et al., 1992). Current integral-based migration algorithms, such as Kirchhoff migrations, are known to be very effective as the ray tracing could be easily adapted to the anisotropic Eikonal equation (Han and Xu, 2012). These algorithms are developed applicable for a single mode of waves such as quasi-P wave, which allows the migration algorithms to be linearized so as to reduce computational cost. However, for wave-equation-based algorithms, an efficient and accurate simulation of a single mode wave, such as quasi-P or quasi-S wave, is essential to the development of successful migration algorithms, such as reverse time migration (RTM) (Baysal et al., 1983; McMechan, 1983; Whitmore, 1983).
Historically, there are different ways for numerical simulation of single mode of waves. First of these is to solve the full elastic wave equations and then to split the wavefields to separate out the quasi-P wave for further analysis. The wavefield separation may be effective for the isotropic case (Sun et al., 2004), but it is not an easy task for heterogeneous anisotropic cases and might be extremely computationally inefficient (Yan and Sava, 2009; Cheng and Fomel, 2013).
Another way to reduce the computational cost while still maintaining the transverse isotropic (TI) anisotropic wave propagation is to apply an acoustic approximation (Alkhalifah, 2000) to the TI equation by setting the shear wave velocity along the symmetry axis to zero. This leads to a scalar fourth-order differential equation. However, Alkhalifah (2000) does not propose any workable numerical technique for obtaining the wave field from this fourth-order differential equation. Zhou et al. (2006) decompose this fourth-order differential equation into a coupled system of 2×2 second-order differential equations. This results in a more computationally efficient scheme than the original elastic equations. However, applications to the TTI media with variable TTI symmetric axis, especially with existence of abrupt changes on the TTI symmetric axis, demonstrate that this system is not numerically stable: the weak instabilities arise and the noises grow linearly in time (Liu et al., 2009, Zhang et al., 2011). These instabilities have been well analyzed and solutions have been proposed (Bakker and Duveneck, 2011; Zhang et al., 2011, Bube et al., 2012). Although the instabilities could be resolved, setting zero the shear wave velocity along the symmetry axis doesn't rule out the existence of the pseudo-shear waves, which are the intrinsic solutions of this couple system. People have tried various ways to attenuate them (Zhang et al., 2009; Guan et al. 2011); or simply ignore the presence of these shear wave artifacts as current industry practices of TTI RTM, with the expectation that they would be canceled out by imaging conditions and migration stack.
Another way is to compute directly the decoupled quasi-acoustic P wave equation. Given constant anisotropic parameters ε and δ, the decoupled equation could be solved with pseudo-spectral algorithm (Etgen and Brandsberg-Dahl, 2009). Therefore, it allows us to take a set of solutions of different constant ε and δ, proceeding with an interpolation scheme, to numerically solve the quasi-P wave with anisotropic parameters that are mildly changing in spatial domain (Chu et al. 2013). For VTI (transverse isotropic with a vertical symmetry axis) media, the interpolation scheme is two dimensional with respect to 249  and δ. But for TTI media, if the symmetrical axis also varies spatially, two additional dimensions are necessary for the interpolation. This pseudo-spectral plus interpolation method may need to compute a huge set of combination of anisotropic parameters and can be very inefficient when it is applied to a very complicated model. When the complexities of the anisotropic model are mild, the decoupled equation may be solved by approximately optimized separable method (Liu et al., 2009), which can be considered as a special case of low rank approximation (Fomel et al., 2012). The computational cost may be lower than pseudo-spectral plus interpolation method, but still very inefficient and inaccurate.